//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : WavesComputer.cpp
//----------------------------------------------------------------------------

// upwind is not implemented correctly so upwind_ok is not defined!

#include <math.h>
#include <stdio.h>
#include <string>
#include <algorithm> // std::max and std::min
#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "defines.h"
#include "Array1D.h"
#include "Array2D.h"
#include "FloatUtils.h"
#include "IncomingWaves.h"
#include "WeakRefl.h"
#include "LowPassFilter.h"
#include "Solver.h"
#include "SolverFactory.h"
#include "WavesComputer.h"


#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#endif


#define SzzEdef() ( U1E * m_dx1twice)
#define SzzWdef() (-U1W * m_dx1twice)
#define SzzNdef() ( U2N * m_dx2twice)
#define SzzSdef() (-U2S * m_dx2twice)

#define SffEdef() ( U1E * m_dx1twice)
#define SffWdef() (-U1W * m_dx1twice)
#define SffNdef() ( U2N * m_dx2twice)
#define SffSdef() (-U2S * m_dx2twice)


// Added by M.A. de Jong
#include <sys/time.h>
#define START_CPU_TIMER(time, iter) \
{ \
    gettimeofday(&start, NULL); \
    for (int t = 0; t < iter; t++) { \

#define STOP__CPU_TIMER(time, iter) \
    } \
    gettimeofday(&stop,  NULL); \
    time = (stop.tv_sec - start.tv_sec) + \
	 1e-6 * (stop.tv_usec - start.tv_usec); \
    time /= iter; \
}

static REAL tcg = 0;

/// constructor
WavesComputer::WavesComputer(int Nx1, ///< number of data points in x1 direction
                             int Nx2, ///< number of data points in x2 direction
                             const Array2D<REAL> &h0, ///< array with depth data
                             const Array2D<REAL> &U1, ///< array with speed data in x1 direction
                             const Array2D<REAL> &U2, ///< array with speed data in x2 direction
                             const Array2D<SHORT> &wet, ///< array with wet/dry data
                             Array2D<REAL>* pzeta, ///< array with waveheight data
                             Array2D<REAL>* pzetam1, ///< array with previous waveheight data
                             Array2D<REAL>* pphi, ///< array with speed potential data
                             Array2D<REAL>* pphim1, ///< array with previous speed potential data
                             Initialization init, ///< perform wave initialisation?
                             REAL time, ///< start time if not zero
                             REAL deltat, ///< last timestep if not zero start time
                             REAL pstart, ///< [s] fade-in time of pulses and incoming waves
                             const char *solvername, ///< name of solver to use
                             REAL dx1, ///< grid size in x1 direction
                             REAL dx2, ///< grid size in x2 direction
                             REAL g, ///< gravity
                             REAL tau_brk, ///< wave breaking parameter
                             int wvin_spec_count, ///< number of specifications for incoming waves
                             const Array1D<REAL> &wvin_spec_ampl, ///< incoming waves amplitudes
                             const Array1D<REAL> &wvin_spec_phase, ///< incoming waves phases
                             const Array1D<REAL> &wvin_spec_freq, ///< incoming waves frequencies
                             const Array1D<REAL> &wvin_spec_dir, ///< incoming waves directions
                             int bcN, ///< boundary condition north
                             int bcE, ///< boundary condition east
                             int bcS, ///< boundary condition south
                             int bcW, ///< boundary condition west
                             int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                             REAL k_wrbc, ///< weak reflection parameter
                             VariationalModel varmod, ///< variational model type
                             REAL kappa_c, ///< wave number
                             REAL psitol, ///< iterative solver tolerance
                             bool short_wave, ///< short wave equations?
                             bool upwind, ///< upwind scheme?
                             REAL eps, ///< ???
                             int M_lp) ///< low pass filter half width
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_dx1(dx1)
, m_dx2(dx2)
, m_dx1square(FloatUtils::square(dx1))
, m_dx2square(FloatUtils::square(dx2))
, m_dx1twice(1.0 / (2 * dx1))
, m_dx2twice(1.0 / (2 * dx2))
, m_g(g)
, m_tau_brk(tau_brk)
, m_bcN(bcN)
, m_bcE(bcE)
, m_bcS(bcS)
, m_bcW(bcW)
, m_short_wave(short_wave)
, m_eps(eps)
, m_upwind(upwind)
, m_lp(M_lp)
, m_h0(h0)
, m_U1(U1)
, m_U2(U2)
, m_wet(wet)
, m_zeta(*pzeta)
, m_zetam1(*pzetam1)
, m_phi(*pphi)
, m_phim1(*pphim1)
, m_pstart(pstart)
, m_time(time)
, m_dt_old(deltat)

, m_pSolver(0)
, m_pi(4.0 * atan(1.0))
, m_dt_new(1)
, m_zeta_mean(0.0)
, m_sum_en_wgt(0.0)
, m_zetam2(COLUMNS, ROWS) // 
, m_phim2(COLUMNS, ROWS) //
, m_psi(COLUMNS, ROWS) //
, m_psim1(COLUMNS, ROWS) //
, m_psim2(COLUMNS, ROWS) //
, m_dt_phi(COLUMNS, ROWS) //
, m_dt_phim1(COLUMNS, ROWS) //
, m_dt_zeta(COLUMNS, ROWS) //
, m_dt_zetam1(COLUMNS, ROWS) //
, m_matS(COLUMNS, ROWS) //
, m_matC(COLUMNS, ROWS) //
, m_matW(COLUMNS, ROWS) //
, m_bpsi(COLUMNS, ROWS) //
, m_alfa_brk(COLUMNS, ROWS) //
, m_zeta_std(COLUMNS, ROWS) //
, m_SffC(COLUMNS, ROWS) // 
, m_SRCp(COLUMNS, ROWS) // 
, m_SRCz(COLUMNS, ROWS) // 
, m_SppC(COLUMNS, ROWS) // 
, m_SzzC(COLUMNS, ROWS) // 
, m_vD0(COLUMNS, ROWS) // 
, m_en_wgt(COLUMNS, ROWS) //
, m_cN(COLUMNS) // cN(0:Nx1+2)
, m_cS(COLUMNS) // cS(0:Nx1+2)
, m_stoN(COLUMNS) // stoN(0:Nx1+2)
, m_stoS(COLUMNS) // stoS(0:Nx1+2)
, m_zeta_in_N_ext(COLUMNS) // zeta_in_N_ext(0:Nx1+2)
, m_zeta_in_S_ext(COLUMNS) // zeta_in_S_ext(0:Nx1+2)
, m_zeta_in_N_cen(COLUMNS) // zeta_in_N_cen(0:Nx1+2)
, m_zeta_in_S_cen(COLUMNS) // zeta_in_S_cen(0:Nx1+2)
, m_zeta_in_N_int(COLUMNS) // zeta_in_N_int(0:Nx1+2)
, m_zeta_in_S_int(COLUMNS) // zeta_in_S_int(0:Nx1+2)
, m_phi_in_N_ext(COLUMNS) // phi_in_N_ext(0:Nx1+2)
, m_phi_in_S_ext(COLUMNS) // phi_in_S_ext(0:Nx1+2)
, m_phi_in_N_cen(COLUMNS) // phi_in_N_cen(0:Nx1+2)
, m_phi_in_S_cen(COLUMNS) // phi_in_S_cen(0:Nx1+2)
, m_phi_in_N_int(COLUMNS) // phi_in_N_int(0:Nx1+2)
, m_phi_in_S_int(COLUMNS) // phi_in_S_int(0:Nx1+2)
, m_psi_in_N_ext(COLUMNS) // psi_in_N_ext(0:Nx1+2)
, m_psi_in_S_ext(COLUMNS) // psi_in_S_ext(0:Nx1+2)
, m_psi_in_N_cen(COLUMNS) // psi_in_N_cen(0:Nx1+2)
, m_psi_in_S_cen(COLUMNS) // psi_in_S_cen(0:Nx1+2)
, m_psi_in_N_int(COLUMNS) // psi_in_N_int(0:Nx1+2)
, m_psi_in_S_int(COLUMNS) // psi_in_S_int(0:Nx1+2)
, m_cE(ROWS) // cE(0:Nx2+2)
, m_cW(ROWS) // cW(0:Nx2+2)
, m_stoE(ROWS) // stoE(0:Nx2+2)
, m_stoW(ROWS) // stoW(0:Nx2+2)
, m_zeta_in_E_ext(ROWS) // zeta_in_E_ext(0:Nx2+2)
, m_zeta_in_W_ext(ROWS) // zeta_in_W_ext(0:Nx2+2)
, m_zeta_in_E_cen(ROWS) // zeta_in_E_cen(0:Nx2+2)
, m_zeta_in_W_cen(ROWS) // zeta_in_W_cen(0:Nx2+2)
, m_zeta_in_E_int(ROWS) // zeta_in_E_int(0:Nx2+2)
, m_zeta_in_W_int(ROWS) // zeta_in_W_int(0:Nx2+2)
, m_phi_in_E_ext(ROWS) // phi_in_E_ext(0:Nx2+2)
, m_phi_in_W_ext(ROWS) // phi_in_W_ext(0:Nx2+2)
, m_phi_in_E_cen(ROWS) // phi_in_E_cen(0:Nx2+2)
, m_phi_in_W_cen(ROWS) // phi_in_W_cen(0:Nx2+2)
, m_phi_in_E_int(ROWS) // phi_in_E_int(0:Nx2+2)
, m_phi_in_W_int(ROWS) // phi_in_W_int(0:Nx2+2)
, m_psi_in_E_ext(ROWS) // psi_in_E_ext(0:Nx2+2)
, m_psi_in_W_ext(ROWS) // psi_in_W_ext(0:Nx2+2)
, m_psi_in_E_cen(ROWS) // psi_in_E_cen(0:Nx2+2)
, m_psi_in_W_cen(ROWS) // psi_in_W_cen(0:Nx2+2)
, m_psi_in_E_int(ROWS) // psi_in_E_int(0:Nx2+2)
, m_psi_in_W_int(ROWS) // psi_in_W_int(0:Nx2+2)
, m_phiN(COLUMNS, 3) // phiN(0:Nx1+2,0:2)
, m_psiN(COLUMNS, 3) // psiN(0:Nx1+2,0:2)
, m_phiS(COLUMNS, 3) // phiS(0:Nx1+2,0:2)
, m_psiS(COLUMNS, 3) // psiS(0:Nx1+2,0:2)
, m_QphiN(COLUMNS, 3) // m_QphiN(0:Nx1+2,0:2)
, m_QphiS(COLUMNS, 3) // m_QphiS(0:Nx1+2,0:2)
, m_QpsiN(COLUMNS, 3) // m_QpsiN(0:Nx1+2,0:2)
, m_QpsiS(COLUMNS, 3) // m_QpsiS(0:Nx1+2,0:2)
, m_QzetaN(COLUMNS, 3) // m_QzetaN(0:Nx1+2,0:2)
, m_QzetaS(COLUMNS, 3) // m_QzetaS(0:Nx1+2,0:2)
, m_zetaN(COLUMNS, 3) // zetaN(0:Nx1+2,0:2)
, m_zetaS(COLUMNS, 3) // zetaS(0:Nx1+2,0:2)
, m_phiE(ROWS, 3) // phiE(0:Nx2+2,0:2)
, m_psiE(ROWS, 3) // psiE(0:Nx2+2,0:2)
, m_phiW(ROWS, 3) // phiW(0:Nx2+2,0:2)
, m_psiW(ROWS, 3) // psiW(0:Nx2+2,0:2)
, m_QphiE(ROWS, 3) // m_QphiE(0:Nx2+2,0:2)
, m_QpsiE(ROWS, 3) // m_QpsiE(0:Nx2+2,0:2)
, m_QphiW(ROWS, 3) // m_QphiW(0:Nx2+2,0:2)
, m_QpsiW(ROWS, 3) // m_QpsiW(0:Nx2+2,0:2)
, m_QzetaE(ROWS, 3) // m_QzetaE(0:Nx2+2,0:2)
, m_QzetaW(ROWS, 3) // m_QzetaW(0:Nx2+2,0:2)
, m_zetaE(ROWS, 3) // zetaE(0:Nx2+2,0:2)
, m_zetaW(ROWS, 3) // zetaW(0:Nx2+2,0:2)
{

    initialize(solvername,
               somfeld,
               k_wrbc,
               varmod,
               kappa_c,
               psitol,
               wvin_spec_count,
               wvin_spec_ampl,
               wvin_spec_phase,
               wvin_spec_freq,
               wvin_spec_dir,
               init);

}

/// destructor
WavesComputer::~WavesComputer()
{
    delete m_pSolver;
    delete m_pIncomingWaves;
    delete m_pWrBc;
    delete m_pLpFil;
}

/// initialization
void WavesComputer::initialize(const char *solvername, ///< solver to use
                               int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                               REAL k_wrbc, ///< weak reflections
                               VariationalModel varmod, ///< variational model type
                               REAL kappa_c, ///< wave number
                               REAL psitol, ///< iterative solver tolerance
                               int wvin_spec_count, ///< number of wave parameters
                               const Array1D<REAL> &wvin_spec_ampl, ///< wave amplitudes
                               const Array1D<REAL> &wvin_spec_phase, ///< waves phases
                               const Array1D<REAL> &wvin_spec_freq, ///< wave frequencies
                               const Array1D<REAL> &wvin_spec_dir, ///< wave directions
                               Initialization init) ///< initialize waves?
{
    /* limit kappa * h_0 to not too small values */
    /* kappa = max( kappa, pi/2 / h0 ) ; */
    Array2D<REAL> kappa(m_zeta.dim1(), m_zeta.dim2());
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            if (m_h0[i1][i2] > 0.0)
            {
                kappa[i1][i2] = std::max(kappa_c, m_pi / 2 / m_h0[i1][i2]);
            }
            else
            {
                kappa[i1][i2] = 1.0;
            }
        }
    }

    /* wave breaking related variables are initialized */
    m_alfa_brk = 0.0;
    m_zeta_std = 0.0;

    /* construct the system matrices for the evolution of the flow equations */
    Array2D<REAL> vN0(m_zeta.dim1(), m_zeta.dim2());
    Array2D<REAL> vM0(m_zeta.dim1(), m_zeta.dim2());

    /* for the moment, we take a regular grid, constant depth and mean velocity field */
    switch (varmod)
    {
        case parabolicVM:
            // NOTE ASR: this leads to very strange results, in the original Fortran code too!
            for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
            {
                for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
                {
                    m_vD0[i1][i2] = 1.0 / 3.0 * m_h0[i1][i2];
                    vM0[i1][i2] = 1.0 / 3.0 * m_h0[i1][i2];
                    vN0[i1][i2] = 2.0 / 15.0 * FloatUtils::third(m_h0[i1][i2]);
                }
            }
            break;

        case coshypVM:
            for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
            {
                for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
                {
                    if (m_h0[i1][i2] > 0.0)
                    {
                        REAL vS0 = sinh(kappa[i1][i2] * m_h0[i1][i2]);
                        REAL vC0 = cosh(kappa[i1][i2] * m_h0[i1][i2]);
                        m_vD0[i1][i2] = vC0 - vS0 / (kappa[i1][i2] * m_h0[i1][i2]);
                        vM0[i1][i2] = 0.5 * kappa[i1][i2] *
                            (vS0 * vC0 - kappa[i1][i2] * m_h0[i1][i2]);
                        vN0[i1][i2] = 0.5 / kappa[i1][i2] *
                            (2.0 * kappa[i1][i2] * m_h0[i1][i2] * FloatUtils::square(vC0)
                             + kappa[i1][i2] * m_h0[i1][i2]
                             - 3.0 * vS0 * vC0);
                    }
                    else
                    {
                        m_vD0[i1][i2] = 0.0;
                        vM0[i1][i2] = 0.0;
                        vN0[i1][i2] = 0.0;
                    }
                }
            }
            break;

        default:
            ASSERT(false, "Unknown variational model");
    }

    /* initialise system matrices */
    /* row and column 0 and m_Nx1+1, resp. 0 and m_Nx2+1, are ghost points */
    Array2D<REAL> matN(m_zeta.dim1(), m_zeta.dim2());
    Array2D<REAL> matE(m_zeta.dim1(), m_zeta.dim2());
    matE = 0.0;
    m_matW = 0.0;
    matN = 0.0;
    m_matS = 0.0;
    m_matC = 1.0;

    m_SzzC = 1.0;
    m_SffC = 1.0;
    m_SppC = 1.0;

    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            REAL U1E = 0.5 * (m_U1[i1    ][i2    ] + m_U1[i1 + 1][i2    ]);
            REAL U1W = 0.5 * (m_U1[i1 - 1][i2    ] + m_U1[i1    ][i2    ]);
            REAL U2N = 0.5 * (m_U2[i1    ][i2    ] + m_U2[i1    ][i2 + 1]);
            REAL U2S = 0.5 * (m_U2[i1    ][i2 - 1] + m_U2[i1    ][i2    ]);

            REAL SzzE = SzzEdef();
            REAL SzzW = SzzWdef();
            REAL SzzN = SzzNdef();
            REAL SzzS = SzzSdef();
            REAL SffE = SffEdef();
            REAL SffW = SffWdef();
            REAL SffN = SffNdef();
            REAL SffS = SffSdef();
            m_SzzC[i1][i2] = -(SzzE + SzzW + SzzN + SzzS);
            m_SffC[i1][i2] = -(SffE + SffW + SffN + SffS);

            REAL mult1 = 1.0 / m_dx1square;
            REAL mult2 = 1.0 / m_dx2square;
            REAL SppE = -0.5 * (vN0[i1][i2] + vN0[i1 + 1][i2 ]) * mult1;
            REAL SppW = -0.5 * (vN0[i1][i2] + vN0[i1 - 1][i2 ]) * mult1;
            REAL SppN = -0.5 * (vN0[i1][i2] + vN0[i1][ i2 + 1]) * mult2;
            REAL SppS = -0.5 * (vN0[i1][i2] + vN0[i1][ i2 - 1]) * mult2;
            m_SppC[i1][i2] = vM0[i1][i2] - (SppE + SppW + SppN + SppS);

            /* Correct for dry nodes */
            if (m_SppC[i1][i2] <= 0.0)
            {
                SppE = 0.0;
                SppW = 0.0;
                SppN = 0.0;
                SppS = 0.0;
                m_SppC[i1][i2] = 1.0;
            }

            /* field equations */
            matE  [i1][i2] = SppE;
            m_matW[i1][i2] = SppW;
            matN  [i1][i2] = SppN;
            m_matS[i1][i2] = SppS;
            m_matC[i1][i2] = m_SppC[i1][i2];
        }
    }

    /* Compute the source terms for a variable current field */
    m_SRCz = 0.0;
    m_SRCp = 0.0;

    REAL SRC_flux_z;
    REAL SRC_flux_p;
    /* first the fluxes in the x1-direction */
    for (int i1 = 0; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            REAL sU1 = 0.0;
            SRC_flux_z = 0.5 * (m_h0[i1][i2] + m_h0[i1 + 1][i2])
                * 0.5 * ((1.0 + sU1) * m_U1[i1][i2] + (1.0 - sU1) * m_U1[i1 + 1][i2]) / m_dx1;
            SRC_flux_p = 0.5 * (m_vD0[i1][i2] + m_vD0[i1 + 1][i2]) * SRC_flux_z;
            m_SRCz[i1][ i2] = m_SRCz[i1][ i2] + SRC_flux_z;
            m_SRCz[i1 + 1][i2] = m_SRCz[i1 + 1][i2] - SRC_flux_z;
            m_SRCp[i1][ i2] = m_SRCp[i1][ i2] + SRC_flux_p;
            m_SRCp[i1 + 1][i2] = m_SRCp[i1 + 1][i2] - SRC_flux_p;
        }
    }

    /* then the fluxes in the x2-direction */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2; ++i2)
        {
            REAL sU2 = 0.0;
            SRC_flux_z = 0.5 * (m_h0[i1][i2] + m_h0[i1][i2 + 1])
                * 0.5 * ((1.0 + sU2) * m_U2[i1][i2] + (1.0 - sU2) * m_U2[i1][i2 + 1]) / m_dx2;
            SRC_flux_p = 0.5 * (m_vD0[i1][i2] + m_vD0[i1][i2 + 1]) * SRC_flux_z;
            m_SRCz[i1][i2 ] = m_SRCz[i1][i2 ] + SRC_flux_z;
            m_SRCz[i1][i2 + 1] = m_SRCz[i1][i2 + 1] - SRC_flux_z;
            m_SRCp[i1][i2 ] = m_SRCp[i1][i2 ] + SRC_flux_p;
            m_SRCp[i1][i2 + 1] = m_SRCp[i1][i2 + 1] - SRC_flux_p;
        }
    }

    // construct the boundary arrays for the incoming waves
    m_pIncomingWaves = new IncomingWaves(m_Nx1, m_Nx2, m_dx1, m_dx2, m_h0, m_U1, m_U2, m_wet,
                                         kappa, m_vD0, vM0, vN0, wvin_spec_count,
                                         wvin_spec_ampl, wvin_spec_phase,
                                         wvin_spec_freq, wvin_spec_dir,
                                         m_g, varmod);

    /* apply boundary conditions */

    /* north */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        if (m_bcN == 0) /* reflective Neumann boundary */
        {
            m_stoN[i1] = matN[i1][m_Nx2];
            m_matC[i1][m_Nx2] = m_matC[i1][m_Nx2] + matN[i1][m_Nx2];
            matN[i1][m_Nx2] = 0.0;
        }
        else if (m_bcN == 1) /* weakly-reflective boundary */
        {
            m_stoN[i1] = matN[i1][m_Nx2];
            matN[i1][m_Nx2] = 0.0;
            m_cN[i1] = std::min(sqrt(m_g * m_h0[i1][m_Nx2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
        m_cN[i1] = std::min(sqrt(m_g * m_h0[i1][m_Nx2]), sqrt(m_g / k_wrbc));
    }
    m_zetaN = 0.0;
    m_phiN = 0.0;
    m_psiN = 0.0;
    m_QzetaN = 0.0;
    m_QphiN = 0.0;
    m_QpsiN = 0.0;

    /* east */
    for (int i2 = 1; i2 <= m_Nx2; ++i2)
    {
        if (m_bcE == 0) /* reflective Neumann boundary */
        {
            m_stoE[i2] = matE[m_Nx1][i2];
            m_matC[m_Nx1][i2] = m_matC[m_Nx1][i2] + matE[m_Nx1][i2];
            matE[m_Nx1][i2] = 0.0;
        }
        else if (m_bcE == 1) /* weakly-reflective boundary */
        {
            m_stoE[i2] = matE[m_Nx1][i2];
            matE[m_Nx1][i2] = 0.0;
            m_cE[i2] = std::min(sqrt(m_g * m_h0[m_Nx1][i2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaE = 0.0;
    m_phiE = 0.0;
    m_psiE = 0.0;
    m_QzetaE = 0.0;
    m_QphiE = 0.0;
    m_QpsiE = 0.0;

    /* south */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        if (m_bcS == 0) /* reflective Neumann boundary */
        {
            m_stoS[i1] = m_matS[i1][1];
            m_matC[i1][1] = m_matC[i1][1] + m_matS[i1][1];
            m_matS[i1][1] = 0.0;
        }
        else if (m_bcS == 1) /* weakly-reflective boundary */
        {
            m_stoS[i1] = m_matS[i1][1];
            m_matS[i1][1] = 0.0;
            m_cS[i1] = std::min(sqrt(m_g * m_h0[i1][1]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaS = 0.0;
    m_phiS = 0.0;
    m_psiS = 0.0;
    m_QzetaS = 0.0;
    m_QphiS = 0.0;
    m_QpsiS = 0.0;

    /* west */
    for (int i2 = 1; i2 <= m_Nx2; ++i2)
    {
        if (m_bcW == 0) /* reflective Neumann boundary */
        {
            m_stoW[i2] = m_matW[1][i2];
            m_matC[1][i2] = m_matC[1][i2] + m_matW[1][i2];
            m_matW[1][i2] = 0.0;
        }
        else if (m_bcW == 1) /* weakly-reflective boundary */
        {
            m_stoW[i2] = m_matW[1][i2];
            m_matW[1][i2] = 0.0;
            m_cW[i2] = std::min(sqrt(m_g * m_h0[1][i2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaW = 0.0;
    m_phiW = 0.0;
    m_psiW = 0.0;
    m_QzetaW = 0.0;
    m_QphiW = 0.0;
    m_QpsiW = 0.0;


    /* compute the system m_matrices and preconditioner */
    m_pSolver = SolverFactory::create(solvername, &m_matC, &m_matS, &m_matW, m_Nx1, m_Nx2, psitol);

    /* weight factor for energy integrations */
    m_en_wgt = 1.0;
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        m_en_wgt[0] [i2] = 0.0;
        m_en_wgt[m_Nx1 + 1][i2] = 0.0;
        m_en_wgt[1] [i2] = 0.5;
        m_en_wgt[m_Nx1] [i2] = 0.5;
    }
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        m_en_wgt[i1][0 ] = 0.0;
        m_en_wgt[i1][m_Nx2 + 1] = 0.0;
        m_en_wgt[i1][1 ] = 0.5;
        m_en_wgt[i1][m_Nx2 ] = 0.5;
    }
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
        {
            m_sum_en_wgt += m_en_wgt[i1][i2];
        }
    }

    m_pLpFil = new LowPassFilter(m_Nx1, m_Nx2, m_lp, m_wet);
    m_pWrBc = new WeakRefl(somfeld, m_Nx1, m_Nx2);

    if (init != noInit)
    {
        if (init == spectrumInit)
        {
            m_pIncomingWaves->startup(m_time - m_dt_old, MAXDEPTH, &m_zetam1, &m_phim1);
            m_pIncomingWaves->startup(m_time, MAXDEPTH, &m_zeta, &m_phi);
        }
        startup();
    }
    else // everything is at rest
    {
        m_time = 0.0;
        m_zeta = 0.0;
        m_zetam1 = 0.0;
        m_phi = 0.0;
        m_phim1 = 0.0;

        m_zetam2 = 0.0;
        m_phim2 = 0.0;
        m_psi = 0.0;
        m_psim1 = 0.0;
        m_psim2 = 0.0;
        m_dt_zetam1 = 0.0;
        m_dt_phim1 = 0.0;
        m_dt_zeta = 0.0;
        m_dt_phi = 0.0;
    }
}


/// compute one simulation step based upon the previous one
void WavesComputer::compute(REAL dt, ///< delta T
                            const Array2D<REAL> &pressurePulse, ///< array with pulse pattern
                            double *pCgTime) ///< time used by solver
{
    ASSERT(pCgTime, "");

    // Added by M.A. de Jong

    m_dt_old = m_dt_new; // old time step
    m_dt_new = dt;
    m_time += m_dt_new; // present time

    REAL zeta_mean = 0.0;
    REAL pslow = pSlow(m_time, m_pstart);

    { //########################################################################

        // shift old values of zeta and phi, and their time derivatives

        shiftZetaPhi();
        // modifies m_zetam2, m_zetam1, m_phim2, m_phim1, m_dt_zetam1, m_dt_phim1 arrays

        // compute the incoming wave conditions along the boundaries
        m_pIncomingWaves->compute(m_time, pslow,
                                  m_zeta_in_N_ext, m_phi_in_N_ext, m_psi_in_N_ext,
                                  m_zeta_in_N_cen, m_phi_in_N_cen, m_psi_in_N_cen,
                                  m_zeta_in_N_int, m_phi_in_N_int, m_psi_in_N_int,
                                  m_zeta_in_E_ext, m_phi_in_E_ext, m_psi_in_E_ext,
                                  m_zeta_in_E_cen, m_phi_in_E_cen, m_psi_in_E_cen,
                                  m_zeta_in_E_int, m_phi_in_E_int, m_psi_in_E_int,
                                  m_zeta_in_S_ext, m_phi_in_S_ext, m_psi_in_S_ext,
                                  m_zeta_in_S_cen, m_phi_in_S_cen, m_psi_in_S_cen,
                                  m_zeta_in_S_int, m_phi_in_S_int, m_psi_in_S_int,
                                  m_zeta_in_W_ext, m_phi_in_W_ext, m_psi_in_W_ext,
                                  m_zeta_in_W_cen, m_phi_in_W_cen, m_psi_in_W_cen,
                                  m_zeta_in_W_int, m_phi_in_W_int, m_psi_in_W_int);
        // modifies all m_* parameters (one dimensional arrays)

        // check for wave breaking if desired
        if (m_tau_brk > 0.0) // m_tau_brk is constant for all threads
        {
            waveBreaking();
            // modifies m_zeta_std and m_alfa_brk
        }


        boundaryZetaPhi();
        // modifies m_zeta and m_phi arrays


        solvePsi(pslow);
        // modifies m_bpsi array, pslow is constant for all threads


        adjustPsi();
        // modifies m_psi array

        

        /* short_wave, normal is: true , set to false for shallow water eq's. */
        { //--------------------------------------------------------------------
            tcg = seconds();
        } //master--------------------------------------------------------------


        if (m_short_wave)
        {
            /* use short wave corrections */
            // bspi is rechterlid, psi is oplossing

            m_pSolver->solve(&m_bpsi, &m_psi);
        }
        else
        {
            /* shallow-water equations */
        }

        
        { //--------------------------------------------------------------------
            double now = seconds();
            *pCgTime = now - tcg;
//printf("%15.12f-%15.12f=%15.12f\n", now, tcg, *pCgTime);
        } //master--------------------------------------------------------------



        // shift old values of psi
        shiftPsi();
        // modifies m_psim2 and m_psim1 arrays


        // apply Neumann boundary conditions to psi
        boundaryPsi();
        // modifies m_psi array


        // compute the time-derivatives of zeta and phi in the domain interior
        timeDerivatives(pslow, pressurePulse);
        // modifies m_dt_zeta and m_dt_phi arrays



        // time stepping of the weakly-reflective boundary conditions
        weaklyReflective();
        // modifies m_zetaN, m_QzetaN, m_phiN, m_QphiN, m_psiN, m_QpsiN and
        // their East, South and West counterpart N*3 dimensional arrays


        // smoothing of dt_zeta and dt_phi if desired, to prevent odd-even decoupling 
        if (m_eps > 0.0) // m_eps is constant for all threads
        {
            smoothZetaPhi();
            // modifies m_dt_zeta and m_dt_phi arrays
        }


        leapFrog();
        // modifies m_zeta and m_phi arrays


        // apply low-pass filter if desired
        if (m_lp > 0) // m_lp is constant for all threads
        {
            m_pLpFil->filter(m_zeta);
            m_pLpFil->filter(m_phi);
            // modifies m_zeta and m_phi arrays
        }


        // compute mean wave height
        int nx1 = m_Nx1 + 1;
        int nx2 = m_Nx2 + 1;
        for (int i1 = 0; i1 <= nx1; ++i1)
        {
            for (int i2 = 0; i2 <= nx2; ++i2)
            {
                zeta_mean += m_en_wgt[i1][i2] * m_zeta[i1][i2];
            }
        } // for ---------------------------------------------------------------

     } // #######################################################################

    m_zeta_mean = zeta_mean / m_sum_en_wgt;
} // compute


/// saves wave field for continued simulation at a later time
bool WavesComputer::save(const char* filename)
{
    FILE *fp = fopen(filename, "w");
    if (!fp)
    {
        return false;
    }
    fprintf(fp, "%.7f\n", getTime());
    fprintf(fp, "%.7f\n", getDeltaTime());
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j <= m_Nx2 + 1; ++j)
        {
            fprintf(fp, "%.20f\n", m_zeta[i][j]);
            fprintf(fp, "%.20f\n", m_zetam1[i][j]);
            fprintf(fp, "%.20f\n", m_phi[i][j]);
            fprintf(fp, "%.20f\n", m_phim1[i][j]);
        }
    }
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            fprintf(fp, "%.20f\n", m_zetaN[i][j]);
            fprintf(fp, "%.20f\n", m_zetaS[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaN[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaS[i][j]);
            fprintf(fp, "%.20f\n", m_phiN[i][j]);
            fprintf(fp, "%.20f\n", m_phiS[i][j]);
            fprintf(fp, "%.20f\n", m_QphiN[i][j]);
            fprintf(fp, "%.20f\n", m_QphiS[i][j]);
            fprintf(fp, "%.20f\n", m_psiN[i][j]);
            fprintf(fp, "%.20f\n", m_psiS[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiN[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiS[i][j]);
        }
    }
    for (int j = 0; j < m_Nx2 + 1; ++j)
    {
        for (int i = 0; i < 3; ++i)
        {
            fprintf(fp, "%.20f\n", m_zetaE[i][j]);
            fprintf(fp, "%.20f\n", m_zetaW[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaE[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaW[i][j]);
            fprintf(fp, "%.20f\n", m_phiE[i][j]);
            fprintf(fp, "%.20f\n", m_phiW[i][j]);
            fprintf(fp, "%.20f\n", m_QphiE[i][j]);
            fprintf(fp, "%.20f\n", m_QphiW[i][j]);
            fprintf(fp, "%.20f\n", m_psiE[i][j]);
            fprintf(fp, "%.20f\n", m_psiW[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiE[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiW[i][j]);
        }
    }
    fclose(fp);
    return true;
}

/// loads wave field from saved file
bool WavesComputer::load(const char* filename)
{
    FILE *fp = fopen(filename, "r");
    if (!fp)
    {
       return false;
    }
#ifdef USEDOUBLE
    const char format[] = "%lf\n";
#else
    const char format[] = "%f\n";
#endif
    REAL time;
    REAL deltat;
    int ignore = 0;
    ignore = fscanf(fp, format, &time);
    ignore = fscanf(fp, format, &deltat);
    ASSERT(fabs(deltat - m_dt_old) < 0.0001, "deltaT in file differs from test parameter");
    m_time = time;
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j <= m_Nx2 + 1; ++j)
        {
            ignore = fscanf(fp, format, &(m_zeta[i][j]));
            ignore = fscanf(fp, format, &(m_zetam1[i][j]));
            ignore = fscanf(fp, format, &(m_phi[i][j]));
            ignore = fscanf(fp, format, &(m_phim1[i][j]));
        }
    }
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            ignore = fscanf(fp, format, &(m_zetaN[i][j]));
            ignore = fscanf(fp, format, &(m_zetaS[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaN[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaS[i][j]));
            ignore = fscanf(fp, format, &(m_phiN[i][j]));
            ignore = fscanf(fp, format, &(m_phiS[i][j]));
            ignore = fscanf(fp, format, &(m_QphiN[i][j]));
            ignore = fscanf(fp, format, &(m_QphiS[i][j]));
            ignore = fscanf(fp, format, &(m_psiN[i][j]));
            ignore = fscanf(fp, format, &(m_psiS[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiN[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiS[i][j]));
        }
    }
    for (int j = 0; j < m_Nx2 + 1; ++j)
    {
        for (int i = 0; i < 3; ++i)
        {
            ignore = fscanf(fp, format, &(m_zetaE[i][j]));
            ignore = fscanf(fp, format, &(m_zetaW[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaE[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaW[i][j]));
            ignore = fscanf(fp, format, &(m_phiE[i][j]));
            ignore = fscanf(fp, format, &(m_phiW[i][j]));
            ignore = fscanf(fp, format, &(m_QphiE[i][j]));
            ignore = fscanf(fp, format, &(m_QphiW[i][j]));
            ignore = fscanf(fp, format, &(m_psiE[i][j]));
            ignore = fscanf(fp, format, &(m_psiW[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiE[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiW[i][j]));
        }
    }
    fclose(fp);
    return true;
}


/// compute initial wave condition
void WavesComputer::startup()
{
    m_dt_new = m_dt_old;

    //  shift old values of zeta and phi
    m_zetam2 = m_zetam1;
    m_phim2 = m_phim1;
    m_zetam1 = m_zeta;
    m_phim1 = m_phi;

    //    for (int w = 0; w < 10; ++w)
    //    {
    //        weaklyReflective();
    //    }
    //    boundaryZetaPhi();

    // solve psi at the previous time step
    solvePsi(1);
    // modifies m_bpsi array

    if (m_short_wave)
    {
        /* use short wave corrections */
        // bspi is rechterlid, psi is oplossing
        m_pSolver->solve(&m_bpsi, &m_psi);
    }
    else
    {
        /* shallow-water equations */
    }

    // shift old values of psi
    m_psim2 = m_psim1;
    m_psim1 = m_psi;

    // compute the time-derivatives of zeta and phi in the domain interior
    Array2D<REAL> dummy(m_zeta.dim1(), m_zeta.dim2());
    dummy = 0.0; // no pulses at all
    timeDerivatives(1, dummy);
    // modifies m_dt_zeta and m_dt_phi arrays

    // compute mean wave height
    REAL zeta_mean = 0.0;
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
        {
            zeta_mean += m_en_wgt[i1][i2] * m_zetam1[i1][i2];
        }
    }
    m_zeta_mean = zeta_mean / m_sum_en_wgt;

} // startup


/// returns current simulation time
REAL WavesComputer::getTime() const
{
    return m_time;
}


/// returns current time step
REAL WavesComputer::getDeltaTime() const
{
    return m_dt_new;
}


/// determines the wave breaking parameters
void WavesComputer::waveBreaking()
{
    // in incoming waves staan sinussen op de gpu
    
    // this is done inside the compute() parallelization:
    // compute wave breaking for the entire field (parallel for)

    /* local mean and standard deviation of the surface */
    /* check for wave breaking */
    REAL eps_brk = m_dt_new / m_tau_brk;
    REAL one_minus_eps_brk = 1.0 - eps_brk;
    REAL brk = eps_brk / m_tau_brk;
    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2) calculations
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            // save some processor time by reducing address calculations
            REAL* palfabrk = &(m_alfa_brk[i1][i2]);

            m_zeta_std[i1][i2] = sqrt(one_minus_eps_brk * FloatUtils::square(m_zeta_std[i1][i2]) +
                                      eps_brk * FloatUtils::square(m_zeta[i1][i2] - m_zeta_mean));
            *palfabrk *= one_minus_eps_brk;
            if (m_wet[i1][i2] > 0 && (m_zeta_std[i1][i2] > 0.18 * (m_h0[i1][i2] + m_zeta_mean)))
            {
                *palfabrk += brk;
            }
        }
    }
    //for-------------------------------------------------------------------
}

/// applies Neumann boundary conditions to zeta and phi
void WavesComputer::boundaryZetaPhi()
{
    // this is done inside the compute() parallelization:
    // compute the boundaries (4 parallel sections)
    // (there is no overlap in the corners)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // m_Nx1 or m_Nx2 calculations per section
    { //====================================================================
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][nx2] = m_zeta[j1][m_Nx2] + m_zeta_in_N_ext[j1] - m_zeta_in_N_cen[j1];
                    m_phi [j1][nx2] = m_phi [j1][m_Nx2] + m_phi_in_N_ext [j1] - m_phi_in_N_cen[j1];
                }
            }
            else if (m_bcN == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][nx2] = m_zetaN[j1][0] + m_zeta_in_N_ext[j1];
                    m_phi [j1][nx2] = m_phiN [j1][0] + m_phi_in_N_ext [j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[nx1][j2] = m_zeta[m_Nx1][j2] + m_zeta_in_E_ext[j2] - m_zeta_in_E_cen[j2];
                    m_phi [nx1][j2] = m_phi [m_Nx1][j2] + m_phi_in_E_ext [j2] - m_phi_in_E_cen[j2];
                }
            }
            else if (m_bcE == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[nx1][j2] = m_zetaE[j2][0] + m_zeta_in_E_ext[j2];
                    m_phi [nx1][j2] = m_phiE [j2][0] + m_phi_in_E_ext [j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][0] = m_zeta[j1][1] + m_zeta_in_S_ext[j1] - m_zeta_in_S_cen[j1];
                    m_phi [j1][0] = m_phi [j1][1] + m_phi_in_S_ext [j1] - m_phi_in_S_cen[j1];
                }
            }
            else if (m_bcS == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][0] = m_zetaS[j1][0] + m_zeta_in_S_ext[j1];
                    m_phi [j1][0] = m_phiS [j1][0] + m_phi_in_S_ext [j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[0][j2] = m_zeta[1][j2] + m_zeta_in_W_ext[j2] - m_zeta_in_W_cen[j2];
                    m_phi [0][j2] = m_phi [1][j2] + m_phi_in_W_ext [j2] - m_phi_in_W_cen[j2];
                }
            }
            else if (m_bcW == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[0][j2] = m_zetaW[j2][0] + m_zeta_in_W_ext[j2];
                    m_phi [0][j2] = m_phiW [j2][0] + m_phi_in_W_ext [j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
    } //====================================================================
}


// solves psi at the previous time step
void WavesComputer::solvePsi(REAL pslow)
{
    // this is done inside the compute() parallelization:
    // 1) compute m_bpsi for the entire field (parallel for)
    // 2) compute the boundaries minus the corners (4 parallel sections)
    // 3) compute the corners to avoid parallel clashes (single)

    // 1) compute m_bpsi for the entire field (parallel for)
    //    m_Nx1*m_Nx2 calculations
    REAL mult1 = 0.25 / m_dx1square;
    REAL mult2 = 0.25 / m_dx2square;
    for (int j1 = 1; j1 <= m_Nx1; ++j1)
    {
        for (int j2 = 1; j2 <= m_Nx2; ++j2)
        {
            if (m_SppC[j1][j2] > 0.0)
            {
                int j1plus = j1 + 1;
                int j1min  = j1 - 1;
                int j2plus = j2 + 1;
                int j2min  = j2 - 1;
                REAL h0  = m_h0[j1][j2];
                REAL vd0 = m_vD0[j1][j2];
                REAL SpfE = (h0 + m_h0[j1plus][j2    ]) * (vd0 + m_vD0[j1plus][j2    ]) * mult1;
                REAL SpfW = (h0 + m_h0[j1min ][j2    ]) * (vd0 + m_vD0[j1min ][j2    ]) * mult1;
                REAL SpfN = (h0 + m_h0[j1    ][j2plus]) * (vd0 + m_vD0[j1    ][j2plus]) * mult2;
                REAL SpfS = (h0 + m_h0[j1    ][j2min ]) * (vd0 + m_vD0[j1    ][j2min ]) * mult2;
                REAL SpfC = -(SpfE + SpfW + SpfN + SpfS);
                m_bpsi[j1][j2] = - SpfE * m_phi[j1plus][j2    ]
                                 - SpfW * m_phi[j1min ][j2    ]
                                 - SpfN * m_phi[j1    ][j2plus]
                                 - SpfS * m_phi[j1    ][j2min ]
                                 - SpfC * m_phi[j1    ][j2    ]
                                 - m_SRCp[j1][j2] * pslow;
            }
            else
            {
                ASSERT(false,"it really gets here");
                m_bpsi[j1][j2] = - m_SRCp[j1][j2] * pslow;
            }
        }
    }
    //for-------------------------------------------------------------------

    // 2) compute the boundaries minus the corners (4 parallel sections)
    // m_Nx1-2 or m_Nx2-2 calculations per section
    int nx1 = m_Nx1 - 1;
    int nx2 = m_Nx2 - 1;
    { //====================================================================
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 0) /* Neumann b.c. */
            {
                for (int j1 = 2; j1 <= nx1; ++j1)
                {
                    m_bpsi[j1][m_Nx2] -= m_stoN[j1] * (m_psi_in_N_ext[j1] - m_psi_in_N_cen[j1]);
                }
            }
            else if (m_bcN == 1) /* Sommerfeld b.c. */
            {
                for (int j1 = 2; j1 <= nx1; ++j1)
                {
                    m_bpsi[j1][m_Nx2] -= m_stoN[j1] * (m_psiN[j1][0] + m_psi_in_N_ext[j1]);
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 0) /* Neumann b.c. */
            {
                for (int j2 = 2; j2 <= nx2; ++j2)
                {
                    m_bpsi[m_Nx1][j2] -= m_stoE[j2] * (m_psi_in_E_ext[j2] - m_psi_in_E_cen[j2]);
                }
            }
            else if (m_bcE == 1) /* Sommerfeld b.c. */
            {
                for (int j2 = 2; j2 <= nx2; ++j2)
                {
                    m_bpsi[m_Nx1][j2] -= m_stoE[j2] * (m_psiE[j2][0] + m_psi_in_E_ext[j2]);
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 0) /* Neumann b.c. */
            {
                for (int j1 = 2; j1 <= nx1; ++j1)
                {
                    m_bpsi[j1][1] -= m_stoS[j1] * (m_psi_in_S_ext[j1] - m_psi_in_S_cen[j1]);
                }
            }
            else if (m_bcS == 1) /* Sommerfeld b.c. */
            {
                for (int j1 = 2; j1 <= nx1; ++j1)
                {
                    m_bpsi[j1][1] -= m_stoS[j1] * (m_psiS[j1][0] + m_psi_in_S_ext[j1]);
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 0) /* Neumann b.c. */
            {
                for (int j2 = 2; j2 <= nx2; ++j2)
                {
                    m_bpsi[1][j2] -= m_stoW[j2] * (m_psi_in_W_ext[j2] - m_psi_in_W_cen[j2]);
                }
            }
            else if (m_bcW == 1) /* Sommerfeld b.c. */
            {
                for (int j2 = 2; j2 <= nx2; ++j2)
                {
                    m_bpsi[1][j2] -= m_stoW[j2] * (m_psiW[j2][0] + m_psi_in_W_ext[j2]);
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
    } //====================================================================

    // 3) compute the corners to avoid parallel clashes (single)
    //    (only 8 calculations)
    { //----------------------------------------------------------------
        /* north */
        if (m_bcN == 0) /* Neumann b.c. */
        {
            m_bpsi[    1][m_Nx2] -= m_stoN[    1] * (m_psi_in_N_ext[ 1] - m_psi_in_N_cen[ 1]);
            m_bpsi[m_Nx1][m_Nx2] -= m_stoN[m_Nx1] * (m_psi_in_N_ext[m_Nx1] - m_psi_in_N_cen[m_Nx1]);
        }
        else if (m_bcN == 1) /* Sommerfeld b.c. */
        {
            m_bpsi[1    ][m_Nx2] -= m_stoN[    1] * (m_psiN[ 1][0] + m_psi_in_N_ext[ 1]);
            m_bpsi[m_Nx1][m_Nx2] -= m_stoN[m_Nx1] * (m_psiN[m_Nx1][0] + m_psi_in_N_ext[m_Nx1]);
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
        /* east */
        if (m_bcE == 0) /* Neumann b.c. */
        {
            m_bpsi[m_Nx1][1    ] -= m_stoE[    1] * (m_psi_in_E_ext[ 1] - m_psi_in_E_cen[ 1]);
            m_bpsi[m_Nx1][m_Nx2] -= m_stoE[m_Nx2] * (m_psi_in_E_ext[m_Nx2] - m_psi_in_E_cen[m_Nx2]);
        }
        else if (m_bcE == 1) /* Sommerfeld b.c. */
        {
            m_bpsi[m_Nx1][1    ] -= m_stoE[    1] * (m_psiE[ 1][0] + m_psi_in_E_ext[ 1]);
            m_bpsi[m_Nx1][m_Nx2] -= m_stoE[m_Nx2] * (m_psiE[m_Nx2][0] + m_psi_in_E_ext[m_Nx2]);
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
        /* south */
        if (m_bcS == 0) /* Neumann b.c. */
        {
            m_bpsi[1    ][1] -= m_stoS[    1] * (m_psi_in_S_ext[ 1] - m_psi_in_S_cen[ 1]);
            m_bpsi[m_Nx1][1] -= m_stoS[m_Nx1] * (m_psi_in_S_ext[m_Nx1] - m_psi_in_S_cen[m_Nx1]);
        }
        else if (m_bcS == 1) /* Sommerfeld b.c. */
        {
            m_bpsi[1    ][1] -= m_stoS[    1] * (m_psiS[ 1][0] + m_psi_in_S_ext[ 1]);
            m_bpsi[m_Nx1][1] -= m_stoS[m_Nx1] * (m_psiS[m_Nx1][0] + m_psi_in_S_ext[m_Nx1]);
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
        /* west */
        if (m_bcW == 0) /* Neumann b.c. */
        {
            m_bpsi[1][1    ] -= m_stoW[    1] * (m_psi_in_W_ext[ 1] - m_psi_in_W_cen[ 1]);
            m_bpsi[1][m_Nx2] -= m_stoW[m_Nx2] * (m_psi_in_W_ext[m_Nx2] - m_psi_in_W_cen[m_Nx2]);
        }
        else if (m_bcW == 1) /* Sommerfeld b.c. */
        {
            m_bpsi[1][1    ] -= m_stoW[    1] * (m_psiW[ 1][0] + m_psi_in_W_ext[ 1]);
            m_bpsi[1][m_Nx2] -= m_stoW[m_Nx2] * (m_psiW[m_Nx2][0] + m_psi_in_W_ext[m_Nx2]);
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    } // master ----------------------------------------------------------
}


/// applies Neumann boundary conditions to psi
void WavesComputer::boundaryPsi()
{
    // this is done inside the compute() parallelization:
    // compute the boundaries (4 parallel sections)
    // (there is no overlap in the corners)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // m_Nx1 or m_Nx2 calculations per section
    { //====================================================================
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][nx2] = m_psi[j1][m_Nx2] + m_psi_in_N_ext[j1] - m_psi_in_N_cen[j1];
                }
            }
            else if (m_bcN == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][nx2] = m_psiN[j1][0] + m_psi_in_N_ext[j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[nx1][j2] = m_psi[m_Nx1][j2] + m_psi_in_E_ext[j2] - m_psi_in_E_cen[j2];
                }
            }
            else if (m_bcE == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[nx1][j2] = m_psiE[j2][0] + m_psi_in_E_ext[j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][0] = m_psi[j1][1] + m_psi_in_S_ext[j1] - m_psi_in_S_cen[j1];
                }
            }
            else if (m_bcS == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][0] = m_psiS[j1][0] + m_psi_in_S_ext[j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[0][j2] = m_psi[1][j2] + m_psi_in_W_ext[j2] - m_psi_in_W_cen[j2];
                }
            }
            else if (m_bcW == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[0][j2] = m_psiW[j2][0] + m_psi_in_W_ext[j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
    } //====================================================================
}


/// time stepping of the weakly-reflective boundary conditions
void WavesComputer::weaklyReflective()
{
    // this is done inside the compute() parallelization:
    // compute the reflections (4 parallel sections)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // m_Nx1 or m_Nx2 calculations per section
    { //====================================================================
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 1)
            {
                Array1D<REAL> zeta_bnd(nx1);
                Array1D<REAL> phi_bnd(nx1);
                Array1D<REAL> psi_bnd(nx1);
                Array1D<REAL> zeta_int(nx1);
                Array1D<REAL> phi_int(nx1);
                Array1D<REAL> psi_int(nx1);
                for (int i1 = 1; i1 <= m_Nx1; ++i1)
                {
                    zeta_bnd[i1] = m_zeta[i1][m_Nx2] - m_zeta_in_N_cen[i1];
                    phi_bnd [i1] = m_phi [i1][m_Nx2] - m_phi_in_N_cen[i1];
                    psi_bnd [i1] = m_psi [i1][m_Nx2] - m_psi_in_N_cen[i1];
                    zeta_int[i1] = m_zeta[i1][m_Nx2 - 1] - m_zeta_in_N_int[i1];
                    phi_int [i1] = m_phi [i1][m_Nx2 - 1] - m_phi_in_N_int[i1];
                    psi_int [i1] = m_psi [i1][m_Nx2 - 1] - m_psi_in_N_int[i1];
                }
                m_pWrBc->compute(m_zetaN, m_QzetaN, zeta_bnd, zeta_int, m_cN, m_Nx1, m_dx2, m_dx1, m_dt_new);
                m_pWrBc->compute(m_phiN, m_QphiN, phi_bnd, phi_int, m_cN, m_Nx1, m_dx2, m_dx1, m_dt_new);
                m_pWrBc->compute(m_psiN, m_QpsiN, psi_bnd, psi_int, m_cN, m_Nx1, m_dx2, m_dx1, m_dt_new);
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 1)
            {
                Array1D<REAL> zeta_bnd(nx2);
                Array1D<REAL> phi_bnd(nx2);
                Array1D<REAL> psi_bnd(nx2);
                Array1D<REAL> zeta_int(nx2);
                Array1D<REAL> phi_int(nx2);
                Array1D<REAL> psi_int(nx2);
                for (int i2 = 1; i2 <= m_Nx2; ++i2)
                {
                    zeta_bnd[i2] = m_zeta[m_Nx1] [i2] - m_zeta_in_E_cen[i2];
                    phi_bnd [i2] = m_phi [m_Nx1] [i2] - m_phi_in_E_cen[i2];
                    psi_bnd [i2] = m_psi [m_Nx1] [i2] - m_psi_in_E_cen[i2];
                    zeta_int[i2] = m_zeta[m_Nx1 - 1][i2] - m_zeta_in_E_int[i2];
                    phi_int [i2] = m_phi [m_Nx1 - 1][i2] - m_phi_in_E_int[i2];
                    psi_int [i2] = m_psi [m_Nx1 - 1][i2] - m_psi_in_E_int[i2];
                }
                m_pWrBc->compute(m_zetaE, m_QzetaE, zeta_bnd, zeta_int, m_cE, m_Nx2, m_dx1, m_dx2, m_dt_new);
                m_pWrBc->compute(m_phiE, m_QphiE, phi_bnd, phi_int, m_cE, m_Nx2, m_dx1, m_dx2, m_dt_new);
                m_pWrBc->compute(m_psiE, m_QpsiE, psi_bnd, psi_int, m_cE, m_Nx2, m_dx1, m_dx2, m_dt_new);
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 1)
            {
                Array1D<REAL> zeta_bnd(nx1);
                Array1D<REAL> phi_bnd(nx1);
                Array1D<REAL> psi_bnd(nx1);
                Array1D<REAL> zeta_int(nx1);
                Array1D<REAL> phi_int(nx1);
                Array1D<REAL> psi_int(nx1);
                for (int i1 = 1; i1 <= m_Nx1; ++i1)
                {
                    zeta_bnd[i1] = m_zeta[i1][1] - m_zeta_in_S_cen[i1];
                    phi_bnd [i1] = m_phi [i1][1] - m_phi_in_S_cen [i1];
                    psi_bnd [i1] = m_psi [i1][1] - m_psi_in_S_cen [i1];
                    zeta_int[i1] = m_zeta[i1][2] - m_zeta_in_S_int[i1];
                    phi_int [i1] = m_phi [i1][2] - m_phi_in_S_int [i1];
                    psi_int [i1] = m_psi [i1][2] - m_psi_in_S_int [i1];
                }
                m_pWrBc->compute(m_zetaS, m_QzetaS, zeta_bnd, zeta_int, m_cS, m_Nx1, m_dx2, m_dx1, m_dt_new);
                m_pWrBc->compute(m_phiS, m_QphiS, phi_bnd, phi_int, m_cS, m_Nx1, m_dx2, m_dx1, m_dt_new);
                m_pWrBc->compute(m_psiS, m_QpsiS, psi_bnd, psi_int, m_cS, m_Nx1, m_dx2, m_dx1, m_dt_new);
            }
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 1)
            {
                Array1D<REAL> zeta_bnd(nx2);
                Array1D<REAL> phi_bnd(nx2);
                Array1D<REAL> psi_bnd(nx2);
                Array1D<REAL> zeta_int(nx2);
                Array1D<REAL> phi_int(nx2);
                Array1D<REAL> psi_int(nx2);
                for (int i2 = 1; i2 <= m_Nx2; ++i2)
                {
                    zeta_bnd[i2] = m_zeta[1][i2] - m_zeta_in_W_cen[i2];
                    phi_bnd [i2] = m_phi [1][i2] - m_phi_in_W_cen[i2];
                    psi_bnd [i2] = m_psi [1][i2] - m_psi_in_W_cen[i2];
                    zeta_int[i2] = m_zeta[2][i2] - m_zeta_in_W_int[i2];
                    phi_int [i2] = m_phi [2][i2] - m_phi_in_W_int[i2];
                    psi_int [i2] = m_psi [2][i2] - m_psi_in_W_int[i2];
                }
                m_pWrBc->compute(m_zetaW, m_QzetaW, zeta_bnd, zeta_int, m_cW, m_Nx2, m_dx1, m_dx2, m_dt_new);
                m_pWrBc->compute(m_phiW, m_QphiW, phi_bnd, phi_int, m_cW, m_Nx2, m_dx1, m_dx2, m_dt_new);
                m_pWrBc->compute(m_psiW, m_QpsiW, psi_bnd, psi_int, m_cW, m_Nx2, m_dx1, m_dx2, m_dt_new);
            }
        } //section---------------------------------------------------------
    } //====================================================================
}


/// time stepping using leap-frog time integration
void WavesComputer::leapFrog()
{
    // this is done inside the compute() parallelization:
    // compute m_zeta and m_phi for the entire field (parallel for)

    REAL beta = m_dt_new / m_dt_old;
    REAL fac1 = FloatUtils::square(beta);
    REAL fac2 = 1.0 - FloatUtils::square(beta);
    REAL fac3 = (beta + 1.0) * m_dt_new;
    REAL fac4 = 0.5 * m_dt_new;

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2)*2 calculations
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            REAL alfabrk = m_alfa_brk[i1][i2];
            m_zeta[i1][i2] = fac1 * m_zetam2[i1][i2] +
                             fac2 * m_zetam1[i1][i2] +
                             fac3 * m_dt_zeta[i1][i2];
            m_phi[i1][i2] = (fac1 * m_phim2[i1][i2] +
                             (fac2 - fac4 * alfabrk) * m_phim1[i1][i2] +
                             fac3 * m_dt_phi[i1][i2]) / (1.0 + fac4 * alfabrk);
        }
    } // for----------------------------------------------------------------
}


/// compute the time-derivatives of zeta and phi in the domain interior
void WavesComputer::timeDerivatives(REAL pslow, const Array2D<REAL> &pressurePulse)
{
    // this is done inside the compute() parallelization:
    // compute m_dt_zeta and m_dt_phi for the entire field (parallel for)

    REAL SfzC = m_g;
    REAL mult1 = 1.0 / m_dx1square;
    REAL mult2 = 1.0 / m_dx2square;

    // m_Nx1*m_Nx2*2 calculations
    for (int j1 = 1; j1 <= m_Nx1; ++j1)
    {
        for (int j2 = 1; j2 <= m_Nx2; ++j2)
        {
            int j1plus = j1 + 1;
            int j1min  = j1 - 1;
            int j2plus = j2 + 1;
            int j2min  = j2 - 1;
            REAL U1E = 0.5 * (m_U1[j1   ][j2   ] + m_U1[j1plus][j2    ]);
            REAL U1W = 0.5 * (m_U1[j1min][j2   ] + m_U1[j1    ][j2    ]);
            REAL U2N = 0.5 * (m_U2[j1   ][j2   ] + m_U2[j1    ][j2plus]);
            REAL U2S = 0.5 * (m_U2[j1   ][j2min] + m_U2[j1    ][j2    ]);

            // phi naar zeta
            REAL h0  = m_h0[j1][j2];
            REAL h0e = mult1 * (h0 + m_h0[j1plus][j2    ]);
            REAL h0w = mult1 * (h0 + m_h0[j1min ][j2    ]);
            REAL h0n = mult2 * (h0 + m_h0[j1    ][j2plus]);
            REAL h0s = mult2 * (h0 + m_h0[j1    ][j2min ]);
            REAL SzfE = 0.5 * h0e;
            REAL SzfW = 0.5 * h0w;
            REAL SzfN = 0.5 * h0n;
            REAL SzfS = 0.5 * h0s;
            REAL SzfC = -(SzfE + SzfW + SzfN + SzfS);

            // psi naar zeta
            REAL vd0 = m_vD0[j1][j2];
            REAL SzpE = -0.25 * h0e * (vd0 + m_vD0[j1plus][j2    ]);
            REAL SzpW = -0.25 * h0w * (vd0 + m_vD0[j1min ][j2    ]);
            REAL SzpN = -0.25 * h0n * (vd0 + m_vD0[j1    ][j2plus]);
            REAL SzpS = -0.25 * h0s * (vd0 + m_vD0[j1    ][j2min ]);
            REAL SzpC = -(SzpE + SzpW + SzpN + SzpS);

            // zeta naar zeta
            REAL SzzE = SzzEdef();
            REAL SzzW = SzzWdef();
            REAL SzzN = SzzNdef();
            REAL SzzS = SzzSdef();

            // phi naar phi
            REAL SffE = SffEdef();
            REAL SffW = SffWdef();
            REAL SffN = SffNdef();
            REAL SffS = SffSdef();

            REAL phie = m_phi[j1plus][j2];
            REAL phiw = m_phi[j1min ][j2];
            REAL phin = m_phi[j1][j2plus];
            REAL phis = m_phi[j1][j2min ];
            REAL phic = m_phi[j1][j2];
            REAL zetac = m_zeta[j1][j2];

            m_dt_zeta[j1][j2] = - SzzE * m_zeta[j1plus][j2]
                                - SzzW * m_zeta[j1min ][j2]
                                - SzzN * m_zeta[j1][j2plus]
                                - SzzS * m_zeta[j1][j2min ]
                                - m_SzzC[j1][j2] * zetac
                                - SzfE * phie
                                - SzfW * phiw
                                - SzfN * phin
                                - SzfS * phis
                                - SzfC * phic
                                - SzpE * m_psi[j1plus][j2]
                                - SzpW * m_psi[j1min ][j2]
                                - SzpN * m_psi[j1][j2plus]
                                - SzpS * m_psi[j1][j2min ]
                                - SzpC * m_psi[j1][j2]
                                - m_SRCz[j1][j2] * pslow;

            m_dt_phi[j1][j2] = - SfzC * zetac
                               - SffE * phie
                               - SffW * phiw
                               - SffN * phin
                               - SffS * phis
                               - m_SffC[j1][j2] * phic
                               - pressurePulse[j1][j2];
        }
    } // for -------------------------------------------------------------------
}


/// shifts old values of zeta and phi, and their time derivatives
void WavesComputer::shiftZetaPhi()
{
    // this is done inside the compute() parallelization:
    // copy lots of values (parallel for)

#ifdef ok
        size_t nr = m_zeta.dim1() * m_zeta.dim2(); // all arrays have the same size
        memcpy(m_zetam2.getElements(), m_zetam1.getElements(), nr);
        memcpy(m_zeta.getElements(), m_zetam1.getElements(), nr);
        memcpy(m_phim2.getElements(), m_phim1.getElements(), nr);
        memcpy(m_phim1.getElements(), m_phi.getElements(), nr);
        memcpy(m_dt_zetam1.getElements(), m_dt_zeta.getElements(), nr);
        memcpy(m_dt_phim1.getElements(), m_dt_phi.getElements(), nr);
        return;
#endif

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2)*6 assignments
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            // save some processor time by reducing address calculations
            REAL *pzetam1 = &(m_zetam1[i1][i2]);
            REAL *pphim1  = &(m_phim1 [i1][i2]);

            m_zetam2   [i1][i2] = *pzetam1;
            *pzetam1            = m_zeta   [i1][i2];
            m_phim2    [i1][i2] = *pphim1;
            *pphim1             = m_phi    [i1][i2];
            m_dt_zetam1[i1][i2] = m_dt_zeta[i1][i2];
            m_dt_phim1 [i1][i2] = m_dt_phi [i1][i2];
        }
    }
    // for------------------------------------------------------------------
}


/// shifts old values of psi
void WavesComputer::shiftPsi()
{
    // this is done inside the compute() parallelization:
    // copy lots of values (parallel for)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2)*2 assignments
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            // save some processor time by reducing address calculations
            REAL *ppsim1  = &(m_psim1[i1][i2]);

            m_psim2[i1][i2] = *ppsim1;
            *ppsim1         = m_psi  [i1][i2];
        }
    }
    // for------------------------------------------------------------------
}


/// psi = 2 * psim1 - psim2 (result in m_psi)
void WavesComputer::adjustPsi()
{
    // this is done inside the compute() parallelization:
    // compute new value of m_psi (parallel for)
    REAL beta  = m_dt_new / m_dt_old;
    REAL beta1 = beta + 1;
    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2) calculations
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            m_psi[i1][i2] = beta1 * m_psim1[i1][i2] - beta * m_psim2[i1][i2];
        }
    }
    // for--------------------------------------------------------------
}


/// smoothing of dt_zeta and dt_phi, to prevent odd-even decoupling
void WavesComputer::smoothZetaPhi()
{
    // this is done inside the compute() parallelization:
    // compute new values of m_dt_zeta and m_dt_phi (parallel for)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;
    double m_eps_minus = 1 - m_eps;

    // (m_Nx1+2)*(m_Nx2+2) calculations
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            m_dt_zeta[i1][i2] = m_eps_minus * m_dt_zeta[i1][i2] + m_eps * m_dt_zetam1[i1][i2];
            m_dt_phi [i1][i2] = m_eps_minus * m_dt_phi [i1][i2] + m_eps * m_dt_phim1 [i1][i2];
        }
    }
    // for--------------------------------------------------------------
}


